
import java.util.Random;

public class IsingRun
{

    public static void main(String[] args)
    {

//        Ising2D ising = new Ising2D(1000,1000, 2.27);
        IsingBits2D ising = new IsingBits2D(18, 6, 2.27);
        System.out.println(ising.getMagnetization());
        for (int iterations = 0; iterations < 1000; iterations++)
        {
            ising.step();
        }
        System.out.println(ising.getMagnetization());
    }
}

class IsingBits2D
{
    
    /*
     *
     *
     *
     *
     */

    /** Neighbors from north:   S >> 8 */
    /** Neighbors from south:   S << 8 */
    /** Neighbors from east:    S << 1 */
    /** Neighbors from west:    S >> 1 */

//S0 is short for "sum equals zero" etc.
//L1 is short for "layer one" etc.
//
////Look at just the first two layers:
//S0= ~(L1 | L2);
//S1= L1 ^ L2;
//S2= L1 & L2;
//
////Now look at the third layer:
//S3 =  L3 & S2;
//S2 = (S2 & ~L3) | (S1 & L3);
//S1 = (S1 & ~L3) | (S0 & L3);
//S0 = S0 & ~L3;
//
////The fourth layer.
//S4 = S3 & L4;
//S3 = (S3 & ~L4) | (S2 & L4);
//S2 = (S2 & ~L4) | (S1 & L4);
//S1 = (S1 & ~L4) | (S0 & L4);
//S0 = S0 & ~L4;
//
//
// 
//

    /*
     * 11111111
     * 10000001
     * 10000001
     * 10000001
     * 10000001
     * 10000001
     * 10000001
     * 11111111
     */
//    private long border = 0xff818181818181ffL;
    private long border = Long.parseLong(
            "11111111" +
            "10000001" +
            "10000001" +
            "10000001" +
            "10000001" +
            "10000001" +
            "10000001" +
            "11111111" +
            "",2);

//    private long borderLeft = 0x0080808080808000L;
//    private long borderLeft = 0x0080808080808000L;

    private long borderRight = 0x0001010101010100L;
    private long borderTop = 0xff00000000000000L;
    private long borderBottom = 0x00000000000000ffL;

    private final int rows;
    private final int cols;

    private final long[][] spins;

    private final int bitsPerElem = 64;

    private final double T;

    private final double[] prob;

    public IsingBits2D(int rows, int cols, double T)
    {
        this.rows = rows;
        this.cols = cols;
        this.T = T;

   
        /* First and last column will hold values from neighbouring cels, same for rows */
        spins = new long[rows/6][cols/6];

        Random rand = new Random();

        for (int i = 0; i < spins.length; i++)
        {
            for (int j = 0; j < spins[i].length; j++)
            {
                spins[i][j] = rand.nextLong(); // all 1111111 in bits
            }
        }

        boundaryConditions();
        for (int i = 0; i < spins.length; i++)
        {
            for (int j = 0; j < spins[i].length; j++)
            {
               printLong(spins[i][j]); 
            }
        }

        prob = new double[5];
        prob[0] = 1;
        prob[1] = 1;
        prob[2] = 1;
        prob[3] = (double)Math.exp(-4/T);
        prob[4] = (double)Math.exp(-8/T);
    }

    private void show()
    {

        for (int i = 0; i < spins.length; i++)
        {
            for (int bitRow = 0; bitRow < 8; bitRow++)
            {
                for (int j = 0; j < spins[i].length; j++)
                {
//                    System.out.print((Long.toBinaryString(spins[i][j])).substring(0+8*bitRow,8+8*bitRow));
//                    System.out.print(" ");
                }
                System.out.println();
            }
            System.out.println("\n");
        }
    }



    private static void printLong(long n)
    {
        String oneLine = binaryLong(n);

        for (int i = 0; i < 8; i++)
        {
            System.out.println(oneLine.substring(0+ 8*i, 7 + 8*i));
        }
        System.out.println();
    }

    private static String binaryLong(long n)
    {
        String num = Long.toBinaryString(n);
        StringBuilder build = new StringBuilder();
        if (num.length() < 64)
        {
            for (int i = num.length(); i < 64; i++)
            {
                build.append("0");
            }
        }
        build.append(num);

        return build.toString();
    }

    private void boundaryConditions()
    {
        for (int i = 1; i < spins.length - 1; i++)
        {
            for (int j = 1; j < spins[i].length - 1; j++)
            {
                // zero current border
                spins[i][j] = spins[i][j] & ~border;
                
                // add left border from the right neighbour
                spins[i][j] = spins[i][j] | ( (borderLeft & spins[i][j+1]) >>> 7);

                spins[i][j] = spins[i][j] | ( (borderRight & spins[i][j-1]) << 7);

                spins[i][j] = spins[i][j] | ( (borderTop & spins[i+1][j]) >>> 56);

                spins[i][j] = spins[i][j] | ( (borderBottom & spins[i-1][j]) << 56);
            }
        }

        for (int i = 0; i < spins.length; i++)
        {
            // zero current border
            spins[i][0] = spins[i][0] & ~border;

            // add left border from the right neighbour
            spins[i][0] = spins[i][0] | ( (borderLeft & spins[i][spins[i].length+1]) >>> 7);

            spins[i][0] = spins[i][0] | ( (borderRight & spins[i][j-1]) << 7);

            spins[i][0] = spins[i][0] | ( (borderTop & spins[i+1][j]) >>> 56);

            spins[i][0] = spins[i][0] | ( (borderBottom & spins[i-1][j]) << 56);
        }
    }

    private int countNeib(int i, int j)
    {
        long spin = spins[i][j];


        long neighTop = spin >>> 8;
        long neighBottom = spin << 8;
        long neighLeft = spin >>> 1;
        long neighRight = spin << 1;

        long sumEquals0;
        long sumEquals1;
        long sumEquals2;
        long sumEquals3;
        long sumEquals4;

        sumEquals0= ~(neighTop | neighBottom);
        sumEquals1= neighTop ^ neighBottom;
        sumEquals2= neighTop & neighBottom;

        //Now look at the third layer:
        sumEquals3 =  neighLeft & sumEquals2;
        sumEquals2 = (sumEquals2 & ~neighLeft) | (sumEquals1 & neighLeft);
        sumEquals1 = (sumEquals1 & ~neighLeft) | (sumEquals0 & neighLeft);
        sumEquals0 = sumEquals0 & ~neighLeft;

        //The fourth layer.
        sumEquals4 = sumEquals3 & neighRight;
        sumEquals3 = (sumEquals3 & ~neighRight) | (sumEquals2 & neighRight);
        sumEquals2 = (sumEquals2 & ~neighRight) | (sumEquals1 & neighRight);
        sumEquals1 = (sumEquals1 & ~neighRight) | (sumEquals0 & neighRight);
        sumEquals0 = sumEquals0 & ~neighRight;


        return 0;
    }

    public void iterate()
    {

    }

    public double getMagnetization()
    {
        int count = 0;
        for (int i = 0; i < spins.length ; i++)
        {
            for (int j = 0; j < spins[i].length ; j++)
            {
                count += Long.bitCount(spins[i][j] & ~border);
            }
        }
        return ((double)count/(rows*cols));
    }

    public void step()
    {

    };
}









class Ising2D
{
    private final int sizeX;
    private final int sizeY;
    private final double temperature;
    private final int[] spins;

    private double[] prob;   //probability of a flip

    private int magnetization;


    public double getMagnetization()
    {
        return (double)magnetization/(sizeX*sizeY);
    }

    public Ising2D(int n, int m, double T)
    {
        sizeX = n;
        sizeY = m;
        temperature = T;
        spins = new int[(sizeY+2)*(sizeX+2)];

        magnetization = sizeX*sizeY;
        for (int i = 0; i < sizeY + 2; i++)
        {
            for (int j = 0; j < sizeX + 2; j++)
            {
                spins[i*(sizeX+2) + j] = 1;
            }
        }
        prob = new double[5];
        prob[0] = 1;
        prob[1] = 1;
        prob[2] = 1;
        prob[3] = (double)Math.exp(-4/temperature);
        prob[4] = (double)Math.exp(-8/temperature);
    }

    private int countNeib(int i, int j)
    {
        int neib = 0;
        int currentSpin = spins[i*(sizeX+2) + j];
        if (spins[i*(sizeX+2) + j-1] == currentSpin) neib++;
        if (spins[i*(sizeX+2) + j+1] == currentSpin) neib++;
        if (spins[(i-1)*(sizeX+2) + j] == currentSpin) neib++;
        if (spins[(i+1)*(sizeX+2) + j] == currentSpin) neib++;
        return neib;
    }


    private final void iterate(int i, int j)
    {
       int neib = countNeib(i, j);
       if (Math.random() < prob[neib])
       {
            spins[i*(sizeX+2) + j] = -spins[i*(sizeX+2)+j];
            magnetization += 2*spins[i*(sizeX+2) + j];
       }
    }

    private void boundaryConditions()
    {
        for (int i = 0; i < sizeX + 1; i++)
        {
            spins[i] = spins[(sizeY-1)*(sizeX+2) + i];
            spins[sizeY*(sizeX+2) + i] = spins[1*(sizeX+2)+ i];
        }
        for (int i = 0; i < sizeX + 1; i++)
        {
            spins[i*(sizeX+2)] = spins[i*(sizeX+2) + sizeX-1];
            spins[i*(sizeX+2) + sizeX] = spins[i*(sizeX+2) + 1];
        }
    }

    final public void step()
    {

        for (int black_white = 0; black_white < 2; black_white++)
        {
            int beginX = black_white;
            for (int i = 1; i < sizeY+1; i++)
            {
                for (int j = 1 + beginX; j < sizeX+1; j = j + 2)
                {
                    iterate(i,j);
                }
                beginX = 1 - beginX;
            }
            boundaryConditions();
        }
    }

}
